library(data.table)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
##
## between, first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggpubr)
## Loading required package: ggplot2
library(Biostrings)
## Loading required package: BiocGenerics
## Warning: package 'BiocGenerics' was built under R version 4.0.5
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, basename, cbind, colnames, dirname, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
## lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:tidyr':
##
## expand
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following objects are masked from 'package:data.table':
##
## first, second
## The following object is masked from 'package:base':
##
## expand.grid
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## The following object is masked from 'package:data.table':
##
## shift
## Loading required package: XVector
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
##
## strsplit
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ tibble 3.1.4 ✔ stringr 1.4.0
## ✔ readr 2.0.1 ✔ forcats 0.5.1
## ✔ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ BiocGenerics::Position() masks ggplot2::Position(), base::Position()
## ✖ dplyr::between() masks data.table::between()
## ✖ Biostrings::collapse() masks IRanges::collapse(), dplyr::collapse()
## ✖ BiocGenerics::combine() masks dplyr::combine()
## ✖ purrr::compact() masks XVector::compact()
## ✖ IRanges::desc() masks dplyr::desc()
## ✖ S4Vectors::expand() masks tidyr::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ S4Vectors::first() masks dplyr::first(), data.table::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks data.table::last()
## ✖ purrr::reduce() masks IRanges::reduce()
## ✖ S4Vectors::rename() masks dplyr::rename()
## ✖ XVector::slice() masks IRanges::slice(), dplyr::slice()
## ✖ purrr::transpose() masks data.table::transpose()
library(foreach)
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
library(stringdist)
##
## Attaching package: 'stringdist'
## The following object is masked from 'package:tidyr':
##
## extract
library(doParallel)
## Loading required package: iterators
Read in Wuhan binding dataset
- Sample of what the data look like.
WUHAN_PEPTIDES_BINDING_SCORES=readRDS("ORIGINAL_PEPTIDES_BINDING_SCORES_V2.rds")
WUHAN_PEPTIDES_BINDING_SCORES %>% head
## # A tibble: 6 × 15
## Peptide Predicted_Bindi… `Thalf(h)` `EL-score` EL_Rank `BA-score` BA_Rank
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 AADLDDFSKQL HLA-A0101,HLA-A… 0.2226,0.… 0.0085,0.… 5.3632… 0.0571,0.… 19.165…
## 2 AAGLEAPFLY HLA-A0101,HLA-A… 0.9656,0.… 0.2728,4e… 0.5874… 0.2922,0.… 0.7888…
## 3 AAGLEAPFLYLY HLA-A0101,HLA-A… 0.9528,0.… 0.0694,1e… 1.6217… 0.2495,0.… 1.0825…
## 4 AAISDYDYY HLA-A0101,HLA-A… 1.0031,0.… 0.3568,4e… 0.4523… 0.3202,0.… 0.6298…
## 5 AAISDYDYYR HLA-A0101,HLA-A… 0.3461,0.… 0.0039,4e… 8.4074… 0.1012,0.… 6.4802…
## 6 AAVYRINW HLA-A0101,HLA-A… 0.3428,0.… 4e-04,0,0… 30.964… 0.0245,0.… 62.656…
## # … with 8 more variables: Ave <chr>, Binder <chr>, nM_41 <chr>, nM <chr>,
## # Aff_Binder <chr>, Length <int>, HydrophobicCount <int>, HydroFraction <dbl>
WUHAN_PEPTIDES_BINDING_SCORES %>% nrow
## [1] 1380
Read in Omicron mutanta binding dataset
- Sample of what the data look like.
OMICRON_PEPTIDES_BINDING_SCORES=readRDS("OMICRON_PEPTIDES_BINDING_SCORES_V2.rds")
OMICRON_PEPTIDES_BINDING_SCORES %>% head
## # A tibble: 6 × 13
## VariantAlignment MT_Predicted_Bind… `MT_Thalf(h)` `MT_EL-score` MT_EL_Rank
## <chr> <chr> <chr> <chr> <chr>
## 1 AEYVNNSY HLA-A0101,HLA-A02… 0.2295,0.1438,… 0.0412,0,1e-0… 2.2565,55.…
## 2 AKSHNITLIW HLA-A0101,HLA-A02… 0.2395,0.1924,… 0.0048,1e-04,… 7.3782,47.…
## 3 ALRITFGGP HLA-A0101,HLA-A02… 0.0843,0.1071,… 0,9e-04,0.004… 73.4615,18…
## 4 APGQTGNIA HLA-A0101,HLA-A02… 0.0906,0.1022,… 5e-04,0,1e-04… 29.3103,52…
## 5 CNDPFLDHK HLA-A0101,HLA-A02… 0.1508,0.093,0… 0.0105,0,1e-0… 4.7481,63.…
## 6 CNDPFLDHKN HLA-A0101,HLA-A02… 0.0855,0.0757,… 2e-04,0,0,0,0… 46,90,95,9…
## # … with 8 more variables: MT_BA-score <chr>, MT_BA_Rank <chr>, MT_Ave <chr>,
## # MT_Binder <chr>, MT_nM_41 <chr>, Predicted_Binding <chr>, MT_nM <chr>,
## # MT_Aff_Binder <chr>
OMICRON_PEPTIDES_BINDING_SCORES %>% nrow
## [1] 81
Read in mutations
MUTATIONS = readRDS("OMICRON_EPITOPE_MUTATIONS.rds")
MUTATIONS %>% mutate(Length = nchar(Peptide))%>% filter(Protein == 'surface glycoprotein') %>% filter(Length== 9)%>% nrow
## [1] 28
#MUTATIONS=MUTATIONS %>% group_by(Peptide, VariantAlignment,start_pos,Mutation) %>% dplyr::summarise(Protein = paste0(Protein, collapse = ","))
join all
MUTATIONS %>% inner_join(WUHAN_PEPTIDES_BINDING_SCORES)%>% nrow
## Joining, by = "Peptide"
## [1] 80
MUTATIONS %>% inner_join(OMICRON_PEPTIDES_BINDING_SCORES) %>% nrow
## Joining, by = "VariantAlignment"
## [1] 80
FULL_MUTATIONS_DT=MUTATIONS %>% inner_join(WUHAN_PEPTIDES_BINDING_SCORES) %>% inner_join(OMICRON_PEPTIDES_BINDING_SCORES)
## Joining, by = "Peptide"
## Joining, by = c("VariantAlignment", "Predicted_Binding")
#FULL_MUTATIONS_DT %>% head %>% DT::datatable()
FULL_MUTATIONS_DT %>% filter(Peptide %in% "WLLWPVTLA")
## # A tibble: 1 × 32
## Peptide VariantAlignment start_pos Protein MutationPos Mutation Hamming
## <chr> <chr> <int> <chr> <chr> <chr> <dbl>
## 1 WLLWPVTLA WLLWPVTLT 55 membrane gl… 63 A63T 1
## # … with 25 more variables: Predicted_Binding <chr>, Thalf(h) <chr>,
## # EL-score <chr>, EL_Rank <chr>, BA-score <chr>, BA_Rank <chr>, Ave <chr>,
## # Binder <chr>, nM_41 <chr>, nM <chr>, Aff_Binder <chr>, Length <int>,
## # HydrophobicCount <int>, HydroFraction <dbl>, MT_Predicted_Binding <chr>,
## # MT_Thalf(h) <chr>, MT_EL-score <chr>, MT_EL_Rank <chr>, MT_BA-score <chr>,
## # MT_BA_Rank <chr>, MT_Ave <chr>, MT_Binder <chr>, MT_nM_41 <chr>,
## # MT_nM <chr>, MT_Aff_Binder <chr>
FULL_MUTATIONS_DT %>% select(Peptide, VariantAlignment) %>% distinct()
## # A tibble: 80 × 2
## Peptide VariantAlignment
## <chr> <chr>
## 1 AEHVNNSY AEYVNNSY
## 2 AKSHNIALIW AKSHNITLIW
## 3 APGQTGKIA APGQTGNIA
## 4 APRITFGGP ALRITFGGP
## 5 CNDPFLGVY CNDPFLDHK
## 6 CNDPFLGVYY CNDPFLDHKN
## 7 DGVYFASTEK DGVYFASIEK
## 8 DSKVGGNYNY DSKVSGNYNY
## 9 EELKKLLEQW EELKKLLEEW
## 10 FCNDPFLGVY FCNDPFLDHK
## # … with 70 more rows
AGRETOPICITY_DT=FULL_MUTATIONS_DT %>% select(Peptide, VariantAlignment, Protein, Predicted_Binding, nM_41, MT_nM_41,Binder,MT_Binder) %>%
separate_rows_(cols = c("Predicted_Binding","nM_41","MT_nM_41","Binder","MT_Binder"),sep=",") %>% mutate(nM_41 = as.numeric(nM_41)) %>% mutate(MT_nM_41 = as.numeric(MT_nM_41))
AGRETOPICITY_DT %>% head
## # A tibble: 6 × 8
## Peptide VariantAlignment Protein Predicted_Binding nM_41 MT_nM_41 Binder
## <chr> <chr> <chr> <chr> <dbl> <dbl> <chr>
## 1 AEHVNNSY AEYVNNSY surface gl… HLA-A0101 17149. 14867. NONBI…
## 2 AEHVNNSY AEYVNNSY surface gl… HLA-A0201 45116. 42972. NONBI…
## 3 AEHVNNSY AEYVNNSY surface gl… HLA-A0202 43770. 41644. NONBI…
## 4 AEHVNNSY AEYVNNSY surface gl… HLA-A0206 42925. 39324. NONBI…
## 5 AEHVNNSY AEYVNNSY surface gl… HLA-A0211 44727. 42694. NONBI…
## 6 AEHVNNSY AEYVNNSY surface gl… HLA-A0301 40402. 37172. NONBI…
## # … with 1 more variable: MT_Binder <chr>
AGRETOPICITY_DT=AGRETOPICITY_DT %>% mutate(Agretopicity = MT_nM_41/nM_41)
SPIKE_AGRETOPICITY=AGRETOPICITY_DT %>% filter(Protein == 'surface glycoprotein') %>% mutate(Length = nchar(Peptide))
SPIKE_AGRETOPICITY=SPIKE_AGRETOPICITY %>% filter(! (Binder == "NONBINDER" & MT_Binder == "NONBINDER"))
SPIKE_AGRETOPICITY %>% select(Peptide,VariantAlignment) %>% distinct() %>% nrow
## [1] 57
SPIKE_AGRETOPICITY %>% group_by(Peptide) %>% dplyr::summarise(n=n())
## # A tibble: 57 × 2
## Peptide n
## <chr> <int>
## 1 AEHVNNSY 3
## 2 APGQTGKIA 2
## 3 CNDPFLGVY 3
## 4 CNDPFLGVYY 3
## 5 DGVYFASTEK 4
## 6 DSKVGGNYNY 3
## 7 FCNDPFLGVY 8
## 8 FCNDPFLGVYY 2
## 9 FGEVFNATRF 6
## 10 FGEVFNATRFASVY 2
## # … with 47 more rows
SPIKE_AGRETOPICITY %>% nrow
## [1] 463
SPIKE_AGRETOPICITY %>% select(Peptide, nM_41, MT_nM_41)%>% distinct() %>% nrow
## [1] 445
ECDF_NM_PLT=SPIKE_AGRETOPICITY %>% select(Peptide, nM_41, MT_nM_41)%>% pivot_longer(cols = c(nM_41,MT_nM_41), names_to = "Dataset",values_to = "Aff")%>%
mutate(Dataset = replace(Dataset, Dataset == "nM_41", "Wuhan"))%>%mutate(Dataset = replace(Dataset, Dataset == "MT_nM_41", "Omicron"))%>%
mutate(Dataset = factor(Dataset, levels = c("Wuhan","Omicron")))%>%
ggplot(aes(x=Aff, color=Dataset))+stat_ecdf(size=1)+theme_pubr(base_size = 22)+grids()+ guides(color = guide_legend(nrow = 2))+ylab("Cumulative Freq. of Peptides")+xlab("Binding Affinity\n(nM)")+geom_vline(xintercept = 500,linetype="dashed",alpha=0.5)
ECDF_NM_LOG10_PLT=SPIKE_AGRETOPICITY %>% select(Peptide, nM_41, MT_nM_41)%>% pivot_longer(cols = c(nM_41,MT_nM_41), names_to = "Dataset",values_to = "Aff")%>%
mutate(Dataset = replace(Dataset, Dataset == "nM_41", "Wuhan"))%>%mutate(Dataset = replace(Dataset, Dataset == "MT_nM_41", "Omicron"))%>%
mutate(Dataset = factor(Dataset, levels = c("Wuhan","Omicron")))%>%
ggplot(aes(x=Aff, color=Dataset))+stat_ecdf(size=1)+theme_pubr(base_size = 22)+grids()+ guides(color = guide_legend(nrow = 2))+ylab("Cumulative Freq. of Peptides")+xlab("Log10 Binding Affinity\n(nM)")+geom_vline(xintercept = 500,linetype="dashed",alpha=0.5)+scale_x_log10()
VIOLIN_NM_PLT=SPIKE_AGRETOPICITY %>% select(Peptide, nM_41, MT_nM_41)%>% pivot_longer(cols = c(nM_41,MT_nM_41), names_to = "Dataset",values_to = "Aff") %>%
mutate(Dataset = replace(Dataset, Dataset == "nM_41", "Wuhan"))%>%mutate(Dataset = replace(Dataset, Dataset == "MT_nM_41", "Omicron"))%>%
mutate(Dataset = factor(Dataset, levels = c("Wuhan","Omicron")))%>%
ggviolin(x="Dataset",y="Aff",color = "Dataset",add="boxplot",trim=TRUE)+scale_y_log10()+theme_pubr(base_size = 22)+stat_compare_means(label="p.signif",label.x.npc = "center")+ylab("Binding Affinity\n(nM)")
BA1: SUPERTYPES ASSOCIATION
AB_SUPERTYPES=fread("HLA_AB_SUPERTYPES.csv") %>% mutate(Allele = gsub("\\*","",Allele))%>% mutate(Allele = paste0("HLA-",Allele))%>% dplyr::rename(HLA_Allele = Allele)
DATA_FOR_ANALYSIS=FULL_MUTATIONS_DT %>% filter(Protein == 'surface glycoprotein')%>% select(Peptide, VariantAlignment, Predicted_Binding, nM_41, MT_nM_41, Binder,MT_Binder) %>% separate_rows_(cols = c("Predicted_Binding","nM_41","MT_nM_41","Binder","MT_Binder"),sep=",")%>% mutate(nM_41 = as.numeric(nM_41), MT_nM_41=as.numeric(MT_nM_41))%>% mutate(Agretopicity = MT_nM_41/nM_41)
HLA_A_B_AGRETOPICITY=DATA_FOR_ANALYSIS %>% dplyr::rename(HLA_Allele=Predicted_Binding)%>% filter(HLA_Allele %in% grep("HLA-A|HLA-B",HLA_Allele,value = T)) %>% inner_join(AB_SUPERTYPES) %>%
filter(! (Binder == "NONBINDER" & MT_Binder == 'NONBINDER'))
## Joining, by = "HLA_Allele"
HLA_CDATASET_AGRETOPICITY=DATA_FOR_ANALYSIS %>% dplyr::rename(HLA_Allele=Predicted_Binding)%>% filter(!HLA_Allele %in% grep("HLA-A|HLA-B",HLA_Allele,value = T))%>% mutate(Supertype = str_extract(HLA_Allele,"HLA-C[0-9][0-9]"))%>%
mutate(Supertype = gsub("HLA\\-","",Supertype))%>%
filter(! (Binder == "NONBINDER" & MT_Binder == 'NONBINDER'))
BA1_SUPERTYPE_PLT = HLA_A_B_AGRETOPICITY %>% rbind(HLA_CDATASET_AGRETOPICITY)%>% dplyr::rename(Wuhan=nM_41, Omicron=MT_nM_41) %>%
ggpaired(cond1 = "Wuhan", cond2 = "Omicron",line.color = "gray")+theme_pubr(base_size = 20)+facet_wrap(~Supertype,ncol=9)+scale_y_log10()+stat_compare_means(label = "p.signif",label.x.npc = "center")+rotate_x_text(angle=45)+ylab("Binding Affinity (nM)")
PRARR plots
DATA_FOR_SUPERTYPE_MUTPOSHLA=FULL_MUTATIONS_DT %>% mutate(Length = nchar(Peptide)) %>% separate_rows_(cols="MutationPos", sep=",") %>% select(Peptide, VariantAlignment, MutationPos, Predicted_Binding, Binder, MT_Binder, nM_41, MT_nM_41,Mutation)%>%
separate_rows_(cols = "Mutation",sep=",")%>%
separate_rows_(cols = c("Predicted_Binding","Binder","MT_Binder","nM_41","MT_nM_41"), sep=",")%>%mutate(MutationPos = as.numeric(MutationPos))%>%
filter(! (Binder == "NONBINDER" & MT_Binder == "NONBINDER"))%>% mutate(nM_41 = as.numeric(nM_41), MT_nM_41=as.numeric(MT_nM_41))%>% mutate(Agretopicity = MT_nM_41/nM_41)%>% dplyr::rename(HLA_Allele=Predicted_Binding)
DATA_FOR_SUPERTYPE_MUTPOSHLA_C=DATA_FOR_SUPERTYPE_MUTPOSHLA%>% filter(!HLA_Allele %in% grep("HLA-A|HLA-B",HLA_Allele,value = T))%>% mutate(Supertype = str_extract(HLA_Allele,"HLA-C[0-9][0-9]"))
DATA_FOR_SUPERTYPE_MUTPOSHLA_A_B=DATA_FOR_SUPERTYPE_MUTPOSHLA%>% filter(HLA_Allele %in% grep("HLA-A|HLA-B",HLA_Allele,value = T)) %>% inner_join(AB_SUPERTYPES)
## Joining, by = "HLA_Allele"
DATA_FOR_SUPERTYPE_MUTPOSHLA_A_B %>% rbind(DATA_FOR_SUPERTYPE_MUTPOSHLA_C)%>% filter(! Mutation == "NA")%>%
filter(Supertype == "B07")%>% mutate(MT_Binder = factor(MT_Binder, levels = c("NONBINDER","BINDER")))%>% arrange(desc(Agretopicity))%>%select(Peptide, VariantAlignment, HLA_Allele, MT_Binder, Mutation, Agretopicity)%>% distinct %>% DT::datatable()
PRARR_PER_MUT_PLT=DATA_FOR_SUPERTYPE_MUTPOSHLA_A_B %>% rbind(DATA_FOR_SUPERTYPE_MUTPOSHLA_C)%>% filter(! Mutation == "NA")%>%
filter(Supertype == "B07")%>% mutate(MT_Binder = factor(MT_Binder, levels = c("NONBINDER","BINDER")))%>% arrange(desc(Agretopicity))%>%select(Peptide, VariantAlignment, HLA_Allele, MT_Binder, Mutation, Agretopicity)%>% distinct %>%
mutate(pMHC = paste0(Peptide,"_", gsub("HLA-","",HLA_Allele)))%>%
ggbarplot(x="pMHC",y="Agretopicity",fill="MT_Binder",position = position_dodge2())+theme_pubr(base_size = 20)+facet_grid(~Mutation,scales="free",space="free")+scale_y_log10()+rotate_x_text()+ylab("Agretopicity (log10)")+
theme(strip.background = element_blank())
FULL_MUTATIONS_DT %>% filter(Peptide %in% grep("PRRA", Peptide, value=T))%>%
select(Peptide, VariantAlignment, Predicted_Binding, nM_41, MT_nM_41, Binder, MT_Binder)%>%
separate_rows_(cols = c("Predicted_Binding", "nM_41", "MT_nM_41", "Binder", "MT_Binder"),sep=",")%>%
filter(! (Binder == "NONBINDER" & MT_Binder == "NONBINDER"))%>%mutate(nM_41=as.numeric(nM_41))%>% mutate(MT_nM_41=as.numeric(MT_nM_41))%>%
mutate(Agretopicity = MT_nM_41/nM_41)%>% mutate(pMHC = paste0(Peptide,"_", gsub("HLA-","",Predicted_Binding)))%>%
mutate(MT_Binder = factor(MT_Binder, levels = c("NONBINDER","BINDER")))%>%
arrange(desc(Agretopicity))%>%
dplyr::rename("Variant Peptide"=VariantAlignment, "HLA Allele"=Predicted_Binding, WT_BA=nM_41, MT_BA=MT_nM_41)%>%DT::datatable()
PRRAR_PMHC_PLT=FULL_MUTATIONS_DT %>% filter(Peptide %in% grep("PRRA", Peptide, value=T))%>%
select(Peptide, VariantAlignment, Predicted_Binding, nM_41, MT_nM_41, Binder, MT_Binder)%>%
separate_rows_(cols = c("Predicted_Binding", "nM_41", "MT_nM_41", "Binder", "MT_Binder"),sep=",")%>%
filter(! (Binder == "NONBINDER" & MT_Binder == "NONBINDER"))%>%mutate(nM_41=as.numeric(nM_41))%>% mutate(MT_nM_41=as.numeric(MT_nM_41))%>%
mutate(Agretopicity = MT_nM_41/nM_41)%>% mutate(pMHC = paste0(Peptide,"_", gsub("HLA-","",Predicted_Binding)))%>% mutate(MT_Binder = factor(MT_Binder, levels = c("NONBINDER","BINDER")))%>%
arrange(desc(Agretopicity))%>%
ggbarplot(x="pMHC",y="Agretopicity",fill="MT_Binder")+theme_pubr(base_size = 20)+rotate_x_text()+scale_y_log10()+ylab("Agretopicity (log10)")
NUM_HLAS_BIND_PEPS_PLT=FULL_MUTATIONS_DT %>% select(Peptide,Binder, Predicted_Binding, MT_Binder)%>% separate_rows_(cols = c("Binder","Predicted_Binding","MT_Binder"),sep=",")%>% filter(! (Binder == "NONBINDER" & MT_Binder == "NONBINDER"))%>%
pivot_longer(cols = c(Binder,MT_Binder))%>% group_by(Peptide, name,value)%>% dplyr::summarise(n=n())%>% filter(! value == "NONBINDER")%>% arrange(desc(name), desc(n))%>%
dplyr::rename(Dataset=name)%>%mutate(Dataset = replace(Dataset, Dataset == "MT_Binder", "Omicron_Binders"))%>%mutate(Dataset = replace(Dataset, Dataset == "Binder", "Wuhan_Binders"))%>%
mutate(Dataset = factor(Dataset, levels = c("Wuhan_Binders","Omicron_Binders")))%>%
ggbarplot(x="Peptide",y="n", fill="Dataset",position = position_dodge())+theme_pubr(base_size = 20)+rotate_x_text()+ylab("# HLAs predicted to bind")
## `summarise()` has grouped output by 'Peptide', 'name'. You can override using the `.groups` argument.
SUBVAR_ECDF_PLT=readRDS("SUBVAR_ECDF_PLT.rds")+ylab("Cumulative Freq. of Peptides")+theme_pubr(base_size = 22)
SUBVAR_VIOLIN_PLT = readRDS("SUBVAR_VIOLIN_PLT.rds")+theme_pubr(base_size = 22)
SUBVAR_SUPERTYPE_BOXPLT = readRDS("SUBVAR_SUPERTYPE_BOXPLT.rds")+theme_pubr(base_size = 22)
BA1_SUPERTYPE_PLT = HLA_A_B_AGRETOPICITY %>% rbind(HLA_CDATASET_AGRETOPICITY)%>% dplyr::rename(Wuhan=nM_41, Omicron=MT_nM_41) %>% filter(Supertype %in% c("A02","A03","B07","C01"))%>%
ggpaired(cond1 = "Wuhan", cond2 = "Omicron",line.color = "gray")+theme_pubr(base_size = 20)+facet_wrap(~Supertype,ncol=9)+scale_y_log10()+stat_compare_means(label = "p.signif",label.x.npc = "center")+rotate_x_text(angle=45)+ylab("Binding Affinity (nM)")
A_GRID=cowplot::plot_grid(ECDF_NM_PLT, ECDF_NM_LOG10_PLT, VIOLIN_NM_PLT, align="hv",nrow=1)
B_GRID_i=cowplot::plot_grid(SUBVAR_ECDF_PLT,SUBVAR_VIOLIN_PLT, align="hv",nrow=1,axis="blt")
B_GRID_ii=cowplot::plot_grid(SUBVAR_SUPERTYPE_BOXPLT, align="hv",nrow=1,axis="blt")
C_GRID = cowplot::plot_grid(PRARR_PER_MUT_PLT, nrow=1, align="hv", axis="bl")
D_F_GRID = cowplot::plot_grid( PRRAR_PMHC_PLT , nrow=1,align="hv",axis="bl")
E_GRID = cowplot::plot_grid( PRRAR_PMHC_PLT+labs(fill=""),NUM_HLAS_BIND_PEPS_PLT , nrow=1, rel_widths = c(0.20,0.8))
V3
cowplot::plot_grid(A_GRID, B_GRID_i,B_GRID_ii, C_GRID, E_GRID, nrow =5, rel_heights = c(0.8,0.7,1.7,1.1,0.9), align="v", axis="blt")
